rep<-function(c0,c1,c2)
{
  
  m1=max(2,floor(sqrt(c0/sqrt(c1))))
  m2=max(2,floor(sqrt(c0/sqrt(c2))))
  n1<-c()
  n2<-c()
  N1<-c()
  N2<-c()
  for(i in 1:1000){
    ppi1<-rgamma(1000,15,0.8)
    ppi2<-rgamma(1000,20,0.8)
    ctr1<-rgamma(1000,14,0.8)
    ctr2<-rgamma(1000,13,0.8)
    ppid<-ppi2-ppi1
    ctrd<-ctr2-ctr1
    
    xi1<-var(ppi1)+var(ppid)+var(ppi2)
    xi2<-var(ctr1)+var(ctrd)+var(ctr2)
    n1<-c(n1,c0/sqrt(c1)*(sqrt(xi1)/(sqrt(xi1)*sqrt(c1)+sqrt(xi2)*sqrt(c2))))
    n2<-c(n2,c0/sqrt(c2)*(sqrt(xi2)/(sqrt(xi1)*sqrt(c1)+sqrt(xi2)*sqrt(c2))))
    #print(n1)
    #print(n2)
    
    
    v1<-var(ppi1[1:m1])+var(ppi2[1:m1])+var(ppid[1:m1])
    v2<-var(ctr1[1:m2])+var(ctr2[1:m2])+var(ctrd[1:m2])
    
    N1=c(N1,max(m1,floor(c0/sqrt(c1)*v1/(v1*sqrt(c1)+v2*sqrt(c2)))))
    N2=c(N2,floor((c0-c1*N1)/c2))
  }
  print(mean(n1))
  print(mean(n2))
  print(mean(N1))
  print(sd(N1)*sqrt(1/length(N1)))
  print(mean(N2))
  print(sd(N2)*sqrt(1/length(N2)))
  
  fppi1<-ppi1[1:mean(N1)]
  fppi2<-ppi2[1:mean(N1)]
  fctr1<-ctr1[1:mean(N2)]
  fctr2<-ctr2[1:mean(N2)]
  fppid<-fppi2-fppi1
  fctrd<-fctr2-fctr1
  
  va1<-var(fppi1)/mean(N1)+var(fctr1)/mean(N2)
  o1<-pnorm(3/sqrt(va1)-1.96)+pnorm(-2.5/sqrt(va1)-1.96)
  va2<-var(fppid)/mean(N1)
  o2<-pnorm(3/sqrt(va2)-1.64)
  del2<-sqrt(va2)*(qnorm(0.9)+qnorm(0.9))
  va3<-var(fctrd)/mean(N2)
  o3<-pnorm(3/sqrt(va3)-1.96)+pnorm(-2.5/sqrt(va3)-1.96)
  va4<-var(fppi2)/mean(N1)+var(fctr2)/mean(N2)
  o4<-pnorm(3/sqrt(va4)-1.64)
  del4<-sqrt(va4)*(qnorm(0.9)+qnorm(0.9))
  print(o1)
  print(o2)
  print(o3)
  print(o4)
  print(del2)
  print(del4)
  
  
  compute_delta_t <- function(power, alpha, var_D) {
    z_alpha_2 <- qnorm(1 - alpha / 2)
    
    solve_delta <- function(delta_t) {
      term1 <- pnorm(delta_t / sqrt(var_D) - z_alpha_2)
      term2 <- pnorm(-delta_t / sqrt(var_D) - z_alpha_2)
      return(term1 + term2 - power)
    }
    
    result <- uniroot(solve_delta, lower = 0, upper = 10)
    return(result$root)
  }
  del1 <- compute_delta_t(0.9, 0.05, va1)
  del3 <- compute_delta_t(0.9, 0.05, va3)
  print(del1)
  print(del3)
}

rep(1000,10,8)
rep(1500,20,15)
rep(2000,8,5)



#hypothesis testing

#H1
t.test(fppi1,fctr1)

#H2
t.test(fppi2,fppi1,
       alternative = "greater",
       mu = 0, 
       paired = TRUE,   
       var.equal = TRUE,
       conf.level = 0.95)

#H3
t.test(fctr2,fctr1,
       alternative = "greater",
       mu = 0, 
       paired = TRUE,   
       var.equal = TRUE,
       conf.level = 0.95)

#H4
t.test(fppi2,fctr2)
